Analyses of swisscom data

Grid preparation

Ancillary data

Municipality boundaries

st_gg21 <- st_read("data-raw/swisstopo/swissboundaries3d_2021-07_2056_5728/SHAPEFILE_LV95_LN02/swissBOUNDARIES3D_1_3_TLM_HOHEITSGEBIET.shp",
                   as_tibble = TRUE) %>% 
  st_zm(drop = TRUE, what = "ZM") %>% 
  # move to LV03
  st_transform(21781) %>%
  # filter(! BFS_NUMMER %in% c(2391, 5391, 5394)) %>% 
  # exclude lakes
  # filter(OBJEKTART != "Kantonsgebiet") %>% 
  # exclude FL & enclaves
  filter(ICC == "CH") %>% 
  select(BFS_NUMMER, NAME, KANTONSNUM, GEM_TEIL) %>% 
  rename(GMDNR = BFS_NUMMER,
         GMDNAME = NAME,
         KTNR = KANTONSNUM) %>% 
  arrange(GMDNR)

write_rds(st_gg21, "data/swisstopo/st_gg21.Rds")

STATPOP offset

Could be used to define offset for st_make_grid

statpop_20 <- read_delim("data-raw/BfS/ag-b-00.03-vz2020statpop/STATPOP2020.zip", 
                         delim = ";", escape_double = FALSE, trim_ws = TRUE)[] %>% 
  mutate_all(as.integer)

# statpop_20 %>% slice(1:10) %>% View()

write_rds(statpop_20, "data/BfS/statpop_20.Rds")

offset_bfs <- read_rds("data/BfS/statpop_20.Rds") %>% 
  summarise(min_x = min(X_KOORD),
            min_y = min(Y_KOORD)) %>% 
  st_as_sf(coords = c("min_x", "min_y"), 
           crs = 21781,
           remove = FALSE)

swisscom offset

Using two communities that are furthest away in southerly and westerly directions:

Chansy to define x
- Postal code 1284
- SFOS number 6611

Chiasso to define y
- Postal code 6830
- SFOS number 5250

Tile definitions were pulled from API using query_postal_codes_heatmaps_api.py.

Points of grid were defined using lower left corner coordinates.

WGS84 coordinates were transformed to LV03.

chansy_grid <- fromJSON("data/swisscom/chansy_grid.json")
chansy_grid <- flatten(chansy_grid$tiles) %>% 
  as_tibble()

chansy_grid_sf <- chansy_grid %>% 
  st_as_sf(coords = c("ll.x", "ll.y"), 
           crs = 4326,
           remove = FALSE) %>% 
  st_transform(21781) %>% 
  mutate(x = st_coordinates(.)[, 1],
         y = st_coordinates(.)[, 2])

chiasso_grid <- fromJSON("data/swisscom/chiasso_grid.json")
chiasso_grid <- flatten(chiasso_grid$tiles) %>% 
  as_tibble()

chiasso_grid_sf <- chiasso_grid %>% 
  st_as_sf(coords = c("ll.x", "ll.y"), 
           crs = 4326,
           remove = FALSE) %>% 
  st_transform(21781) %>% 
  mutate(x = st_coordinates(.)[, 1],
         y = st_coordinates(.)[, 2])

Example of grid

Coordinates of lower left corner were then obtained by getting minimum values

offset_swisscom <- 
  bind_rows(
    
  chansy_grid_sf %>% 
  st_drop_geometry() %>%
  summarise(min_x = min(x),
            min_y = min(y)) ,

  chiasso_grid_sf %>% 
  st_drop_geometry() %>%
  summarise(min_x = min(x),
            min_y = min(y)) 
  ) %>%
  summarise(min_x = min(min_x),
            min_y = min(min_y)) %>% 
  # needs rounding - transfrom error?
  mutate(min_y = round(min_y)) %>% 
  st_as_sf(coords = c("min_x", "min_y"), 
           crs = 21781,
           remove = FALSE)

Note small difference from BfS derived minimums (in black)!

Grid

Generate whole country

100m grid, with lower left corner defined using STATPOP cells.

Also, adding ID for each cell (based on row number).

country <- st_gg21 %>% 
  st_make_grid(cellsize = 100, 
               offset = c(offset_swisscom$min_x, offset_swisscom$min_y),
               square = TRUE) %>% 
  st_sf() %>%
  st_cast("POLYGON") %>% 
  mutate(ID1 = row_number()) %>% 
  relocate(ID1)

# st_write(country, "data/grid/country.gpkg")
write_rds(country, "data/grid/country.Rds")

Generate cities

Clipped / selected using swisstopo municipality data.

Added second, city specific ID and area.

Function with clipping

# why st_intersection 
# https://stackoverflow.com/questions/62442150/why-use-st-intersection-rather-than-st-intersects

city_grid_clip <- function(municipality){
  
  city <- country %>% 
    st_intersection(st_gg21 %>% 
                      filter(GMDNR == municipality))  %>%
    mutate(ID2 = row_number(),
           area = st_area(.)) %>% 
    relocate(ID1, ID2, area)
  
  return(city)
}

Function without clipping

# solution from
# https://stackoverflow.com/questions/57014381/how-to-filter-an-r-simple-features-collection-using-sf-methods-like-st-intersect

city_grid <- function(municipality){
  
  city <- country %>% 
    filter(st_intersects(geometry, 
                         st_gg21 %>% filter(GMDNR == municipality), 
                         sparse = FALSE)) %>%
    mutate(ID2 = row_number()) %>% 
    relocate(ID1, ID2)
  
  return(city)
}

# or even simpler, without dplyr
# could be generalized with arguments to provide full grid and GMDE shape

city_grid <- function(municipality){
  
  city <- country[st_gg21 %>% filter(GMDNR == municipality), ] %>%
    mutate(ID2 = row_number()) %>% 
    relocate(ID1, ID2)
  return(city)
  
}

Zuirich

# zurich_clip <- city_grid_clip(261)
zurich <- city_grid(261)

Geneva

# geneva_clip <- city_grid_clip(6621)
geneva <- city_grid(6621)

Basel

# basel_clip <- city_grid_clip(2701)
basel <- city_grid(2701)

Lausanne

# lausanne_clip <- city_grid_clip(5586)
lausanne <- city_grid(5586)

Bern

bern_clip <- city_grid_clip(351)
bern <- city_grid(351)

Offset

Using two communities to get min x and y coordinates: Chansy & Chiasso

# chansy_clip <- city_grid_clip(6611)
# st_write(chansy_clip, "data/grid/chansy_clip.gpkg")
# write_rds(chansy_clip, "data/grid/chansy_clip.Rds")

chansy <- city_grid(6611)
# st_write(chansy, "data/grid/chansy.gpkg")
write_rds(chansy, "data/grid/chansy.Rds")

# chiasso_clip <- city_grid_clip(5250)
# st_write(chiasso_clip, "data/grid/chiasso_clip.gpkg")
# write_rds(chiasso_clip, "data/grid/chiasso_clip.Rds")

chiasso <- city_grid(5250)
# st_write(chiasso, "data/grid/chiasso.gpkg")
write_rds(chiasso, "data/grid/chiasso.Rds")

Example

Municipality Bern coverage:

Clipped grid:

Vs. unclipped one:

Note, when using clipped grid, there are few examples of grid cells with very small area.
These could be excluded if no population was found there?

       ID1  ID2            area GMDNR GMDNAME KTNR GEM_TEIL
1  4259841  401 0.1571798 [m^2]   351    Bern    2        0
2  4332927 2876 0.1766709 [m^2]   351    Bern    2        0
3  4423647 5188 0.2123485 [m^2]   351    Bern    2        0
4  4339904 3098 0.2957572 [m^2]   351    Bern    2        0
5  4273711  702 0.7733263 [m^2]   351    Bern    2        0
6  4245805  139 1.0971039 [m^2]   351    Bern    2        0
7  4245786  123 1.1456480 [m^2]   351    Bern    2        0
8  4420099 5082 1.9921709 [m^2]   351    Bern    2        0
9  4291217 1374 2.4771009 [m^2]   351    Bern    2        0
10 4273760  735 3.2587531 [m^2]   351    Bern    2        0

Area estimates

Bern - hourly

Number of tiles for area of Bern

(tiles <- nrow(bern))

[1] 5487

Number of requests needed

(requests <- ceiling(tiles / 100))

[1] 55

Price for complete grid

(price <- requests * 0.01)

[1] 0.55

Price for 24h

(price_24h <- price * 24)

[1] 13.2

Price for 4w hourly

(price_4w <- price_24h * 7 * 4)

[1] 369.6

Price for year hourly

(price_year <- price_24h * 365)

[1] 4818

Price for 4w daily

(price_4w_d <- price * 7 * 4)

[1] 15.4

Price for year daily

(price_year_d <- price * 365)

[1] 200.75

5 cities - hourly